#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>

#include "XDM_Mesh/XDM_mesh_triangulation_3.h"
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include "XDM_Mesh/XDM_mesh_criteria_3.h"


#include <CGAL/make_mesh_3.h>
#include <CGAL/Image_3.h>
#include "XDM_Mesh/XDM_mesh_domain_3.h"
#include "XDM_Mesh/XDM_Data.h"
#include "XDM_Mesh/XDM_make_mesh.h"
#include <fstream>


// Domain
struct K: public CGAL::Exact_predicates_inexact_constructions_kernel {};
// typedef CGAL::Image_3 Image;
//typedef CGAL::Labeled_image_mesh_domain_3<Image,K> Mesh_domain;
//typedef CGAL::TestMesh<K> Mesh_domain;

// // Triangulation
// typedef CGAL::Mesh_triangulation_3<Mesh_domain>::type Tr;
// typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;

// Mesh Criteria


typedef unsigned char XDMIndexT;
using std::vector;
typedef vector< vector< vector< XDMIndexT > > >  DataMatrixT;

typedef CGAL::XDM_test::XDM_Data XDM_Data;
typedef CGAL::XDM_mesh_domain<XDM_Data, K> XDM_Mesh_Domain;

typedef CGAL::XDM_mesh_triangulation_3<XDM_Mesh_Domain>::type XDM_Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<XDM_Tr> XDM_C3t3;

typedef CGAL::Mesh_criteria_3<XDM_Tr> Mesh_criteria;
typedef Mesh_criteria::Facet_criteria    Facet_criteria;
typedef Mesh_criteria::Cell_criteria     Cell_criteria;

void ReadDumbFile( XDM_Data & oData )
{
  using std::ifstream;
  ifstream is;
  
  is.open("Tetra-35Arb3Ran.dx");
  if( !is.good() )
  {
    std::cerr << "File not opened " << std::endl;
    exit(0);
  }
  int nZ = 62;
  int nX = 62;
  int nY = 62;
  float cx = 31;
  float cy = 31;
  float cz = 31;
  float r_sq = 961;
  oData.SetPixelSize( 1, 1, 1 );
  oData.Resize( nX, nY, nZ );
  for( int k = 0; k < nZ; k ++ )
  {
    for( int i = 0; i < nX; i ++ )
    {
      int nTmp;
      for( int j = 0; j < nY; j ++ )
      {
//         is >> std::skipws >> nTmp;
//         oData.Set( j, i, k, nTmp );
        float fDx = i - cx;
        float fDy = j - cy;
        float fDz = k - cz;
        if( fDx * fDx + fDy * fDy + fDz * fDz <= r_sq )
        {
          if( i > cx && j >= cy )
            oData.Set(i, j, k, 1);
          if( i > cx && j < cy )
            oData.Set(i, j, k, 2);
          if( i <= cx && j >= cy )
            oData.Set(i, j, k, 3);
          if( i <= cx && j < cy )
            oData.Set(i, j, k, 4);

        }
      }
    }
  }
  is.close();
}

int main()
{
  XDM_Data oTestData;
  
  ReadDumbFile(  oTestData );
  oTestData.ReadBndPoints( "TestBnd.txt" );
  XDM_Mesh_Domain oTestDomain( oTestData );

//   // Mesh criteria
//  Facet_criteria facet_criteria(30, 6, 4); // angle, size, approximation
  Cell_criteria cell_criteria(4, 5); // radius-edge ratio, size
   
  Facet_criteria facet_criteria(20, 9, 3); // angle, size, approximation
   //Cell_criteria cell_criteria(1, 1); // radius-edge ratio, size
  
   Mesh_criteria criteria(facet_criteria, cell_criteria);

  // Meshing
   XDM_C3t3 c3t3 = CGAL::XDM_make_mesh_3<XDM_C3t3>( oTestDomain, criteria);
   std::cout << "done " << std::endl;
  
  // Output
   std::ofstream medit_file("out.mesh");
   c3t3.output_to_medit(medit_file);

  return 0;
}
